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Abstract 

This paper presents a review of the current state-of-the-art of numerical 
methods for nonhnear Dirac (NLD) equation. Several methods are extend- 
edly proposed for the (l-l-l)-dimensional NLD equation with the scalar and 
vector self-interaction and analyzed in the way of the accuracy and the time 
reversibility as well as the conservation of the discrete charge, energy and 
linear momentum. Those methods are the Crank-Nicolson (CN) schemes, 
the linearized CN schemes, the odd-even hopscotch scheme, the leapfrog 
scheme, a semi-implicit finite difference scheme, and the exponential opera- 
tor splitting (OS) schemes. The nonlinear subproblems resulted from the OS 
schemes are analytically solved by fully exploiting the local conservation laws 
of the NLD equation. The effectiveness of the various numerical methods, 
with special focus on the error growth and the computational cost, is illus- 
trated on two numerical experiments, compared to two high-order accurate 
Runge-Kutta discontinuous Galerkin methods. Theoretical and numerical 
comparisons show that the high-order accurate OS schemes may compete 
well with other numerical schemes discussed here in terms of the accuracy 
and the efficiency. A fourth-order accurate OS scheme is further applied 
to investigating the interaction dynamics of the NLD solitary waves under 
the scalar and vector self-interaction. The results show that the interaction 
dynamics of two NLD solitary waves depend on the exponent power of the 
self-interaction in the NLD equation; collapse happens after collision of two 
equal one-humped NLD solitary waves under the cubic vector self-interaction 
in contrast to no collapse scattering for corresponding quadric case. 
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1. Introduction 



As a relativistic wave equation, the Dirac equation provides naturally 
a description of an electron [ij. Following Dirac's discovery of the linear 
equation of the electron, there appears the fundamental idea of nonlinear de- 
scription of an elementary spin-1/2 particle which makes it possible to take 
into account its self-interaction Heisenberg put forward the idea to 

use a nonlinear Dirac (NLD) equation as a possible basis model for a unified 
field theory jsj. A key feature of the NLD equation is that it allows solitary 
wave solutions or particle-like solutions — the stable localized solutions with 
finite energy and charge j6|. That is, the particles appear as intense local- 
ized regions of field which can be recognized as the basic ingredient in the 
description of extended objects in quantum field theory 0]. Different self- 
interactions give rise to different NLD models mainly including the Thirring 
model j^, the Soler model joj, the Gross-Neveu model [lo[ (equivalent to 



the massless Soler model), and the bag model ll| {i.e. the solitary waves 
with finite compact support), all of which attracted wide interest of physi- 
cists and mathematicians around the 1970s and 1980s, especially on looking 
for the solitary wave solutions and investigating the related physical and 
mathematical properties j6|. 

For the NLD equation in (1+1) dimensions {i.e. one time dimension plus 
one space dimension), several analytical solitary wave solutions are derived 



m 



12 




for the quadric nonlinearity, [lj| for fractional nonlinearity as 



well as [15|, Il6| for general nonlinearity by using explicitly the constraint 
resulting from energy-momentum conservation, and summarized by Mathieu 
In contrast, there are few explicit solutions in (1+3) dimensions except 



17 



for some particular cases shown in [18[ in spite of their existence claimed by 



mathematicians for various situations 19|-|26| (the readers are referred to an 



overview 



271 ] on this topic), and most understanding is based on numerical 



investigations, e.g. 28-3^1. Beyond this, the study of the NLD equation in 
(1+1) dimensions could be very helpful for that in (1+3) dimensions since the 
(l+l)-dimensional NLD equation correspond to the asymptotic form of the 
equation in the physically interesting case of (1+3) dimensions as emphasized 



by Kaus j3l| • That is, some qualitative properties of the NLD solitary waves 



2 



could be similar in such two cases. An interesting topic for the NLD equation 
is the stability issue, which has been the central topic in works spread out 
over several decades in an effort that is still ongoing. Analytical studies of 



the NLD solitary wave stability face serious obstacles |32|-|34 
of computer simulations are contradictory 



while results 



30, 35-37 



The stability analysis 
of the NLD solitary waves is still a very challenging mathematical problem to 
date 16,381. Recent efforts in this direction can be found in 39-441. Another 



rising mathematical interest related to the NLD equation is the analysis of 
global well-posedness, e.g. see 45|, |46|| and references therein. 

In the case of that theoretical methods were not able to give the satis- 
factory results, numerical methods were used for obtaining the solitary wave 
solutions of the NLD equation as well as for investigating the stability. An 
important step in this direction was made by Alvarez and Carreras in 1981 
7|, who simulated the interaction dynamics between the (l+l)-dimensional 



NLD solitary waves of different initial charge for the Soler model [9| by us- 
ing a second-order accurate Crank-Nicholson (CN) scheme j48|. They first 
saw there: charge and energy interchange except for some particular ini- 
tial velocities of the solitary waves; inelastic interaction in binary collisions; 
and oscillating state production from binary collisions. Motivated by their 
work, Shao and Tang revisited this interaction dynamics problem in 2005 49 
by employing a fourth-order accurate Runge-Kutta discontinuous Galerkin 
(RKDG) method . They not only recovered the phenomena mentioned 
by Alvarez and Carreras but also revealed several new ones, e.g. collapse 
in binary and ternary collisions of two-humped NLD solitary waves 



49 



long-lived oscillating state formed with an approximate constant frequency in 
collisions of two standing waves 53]; full repulsion in binary and ternary col- 



lisions of out-of-phase waves [5l[. Their numerical results also inferred that 
the two-humped profile could undermine the stability during the scattering 
of the NLD solitary waves. Note in passing that the two-humped profile was 



first pointed out by Shao and Tang (49| and later gotten noticed by other 



researchers 16 . Besides the often-used CN 48 



50 



and RKDG methods 
there exist many other numerical schemes for solving the (l+l)-dimensional 



NLD equation: split-step spectral schemes [52 



the semi-implicit scheme |54l] |46 



53|, 

55|. multi-symplectic Runge-Kutta methods 
57| etc. The fourth-order accurate RKDG method 



the linearized CN scheme 
Legendre rational spectral methods 
adap tive mesh methods 
is very appropriate 



56 



50 



for investigating the interaction dynamics of the NLD solitary waves due to 
their ability to capture the discontinuous or strong gradients without pro- 
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ducing spurious oscillations, and thus performs better than the second-order 



accurate CN scheme [48|. However, the high cost due to the relatively more 
freedoms used in each cell and the stringent time step constraint reduce its 
practicality in more realistic simulations where realtime and quantitative re- 
sults are required. 

Recently, there has been a remarkable upsurge of the interest in the NLD 
models, as they emerge naturally as practical models in other physical sys- 



tems, such as the gap solitons in nonlinear optics [38|, Bose-Einstein conden- 
sates in honeycomb optical lattices |58| . as well as matter influencing the evo- 
lution of the Universe in cosmology |59|. In view of such new trend, longtime 
stable, efficient, conservative and high-order accurate numerical methods for 
solving the NLD equation are highly desirable. Finite difference methods, 
usually as the first try in practice, enable easy coding and debugging and 
thus are often used by physicists and engineers. However, detailed discussion 
and careful comparison on finite difference solvers for the NLD equation are 
not existed. To this end, the present work as the first step will extendedly 
propose the finite deference schemes for solving the NLD equation with the 
scalar and vector self-interaction. A general and precise comparison among 
them will be presented. However, all of these finite difference methods are 
often of the second order accuracy and thus sustain fast error growth with 
respect to time. To achieve relatively slow error growth, high-order accurate 
numerical methods are required. By exploiting the local conservation laws of 
the NLD equation, we present exponential operator splitting (OS) schemes 
which are time reversible and can conserve the discrete charge. One of the 
high-order accurate OS schemes is afterwards adopted to investigate the in- 
teraction dynamics for the NLD solitary waves under the general scalar and 
vector self-interaction. It should be noted that the experiments carried out 
in the literatures are all limited to the collisions of the NLD solitary waves 
under the quadric scalar self-interaction. Here, the binary collisions of the 
NLD solitary waves under the cubic scalar self-interaction or under the vec- 
tor self-interaction or under the linear combination of the scalar and vector 
self-interactions are all studied for the first time. 

The paper is organized as follows. There is a brief review of the NLD 
equation in Section |2] and the solitary wave solutions are also derived there 
for the general scalar and vector self-interaction. The numerical schemes 
are presented in Section |3] and corresponding numerical analysis is given in 
Section HJ The numerical results are presented with discussion in Section HJ 
The paper is concluded in Section [H] with a few remarks. 
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2. Preliminaries and notations 



This section will introduce the (l+l)-dimensional NLD equation with the 
scalar and vector self-interaction and derive its two solitary wave solutions. 
Throughout the paper, units in which both the speed of light and the reduced 
Planck constant are equal to one will be used. 

2.1. Nonlinear Dirac equation 

Our attention is restricted to the NLD equation in (1 + 1) dimensions 
which can be written in the covariant form 

{iY^^,-m)^ + dLi[^]/d^ = 0, (2.1) 

where * is the spinor with two complex components, := ^'t^'' denotes the 
adjoint spinor, is the complex conjugate transpose of denotes the 

self- interaction Lagrangian, the Greek index /x runs from to 1, the Einstein 
summation convection has been applied, i is the imaginary unit, m is the 
rest mass, 9^ = -rf-^; stands for the covariant derivative, and 7^^, for u = 0, 1, 



are the gamma matrices or the Dirac matrices, chosen as those in [47|, |49| , 





In fact, Eq. (12. ip is the equation of motion for the classical spinorial parti- 
cle with the Lagrangian being a sum of the Dirac Lagrangian and the self- 
interaction Lagrangian, i.e. 

L[^] = ^{iYd^-m)^ + Li[^]. (2.2) 

There exist several different NLD models in the literature, where two im- 
portant models in (H-l) dimensions are the scalar self-interaction of Soler 

3 

Ls[*] := (2.3) 



and the vector self-interaction of Thirring [60 



L^[q,] ■= *7'^**7^*, (2.4) 

where 7^ = Tj^^'^'^ with the Minkowski metric rj^y = diag(l, —1) on spacetime, 
which implies 7^ = (—1)^7'^. 
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This paper will focus on the NLD equation fl2.ll) with a more general 



self-interaction [15|, |61 



Li[*] = s (L^M)'^' + V {L^l^]^-'^/' , (2.5) 

and extendedly propose and compare its numerical methods, where s and v 
are two real numbers and /c is a positive real number. If the spinor * is scaled 
by a constant factor as = ^/a^ with a G C, then the scaled self-interaction 
Lagrangian will be a'^+^Li[*] which shows that the power exponent to a is 
k + 1. In such sense, we call that the self- interaction Lagrangian Li has the 
power exponent /c-|-l 14l-ll6|. Hereafter the gwadnc and cit6zc self-interaction 
will be referred to the case k = 1 and the case k = 2, respectively. 



The self-interaction ( 12. 5 p implies the so-called homogeneity relation 15 

13 



Combining it with the definition of the Lagrangian L[^] and (12. ip gives 

= -kLi[^]. (2.6) 

For the NLD equation (12. ip with (12. 5p . one may still verify the following 
local conservation laws for the current vector and the energy-momentum 
tensor T^^: 

d'^J^ = 0, df^T^, = 0, (2.7) 

where 

For localized solitary waves ^ = {ipi, 'ip2)'^, one may derive a direct corol- 
lary of (12.70 . i.e. the following global conservation laws 5ol . 

Proposition 2.1. Assume that lim \ilji{x,t)\ = and lim \dxipi{x,t)\ < 

\x\—^+oo \x\—>+oo 

+00 hold uniformly for t > andi = 1, 2. The total energy E, the total linear 
momentum P, and the total charge Q, defined respectively by 

Toodx, P{t) := / Toidx, Q(t) := / Jodx, (2.8) 



satisfy 
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The properties fl2.6p and fl2.7p will be also exploited to find the solitary 
wave solutions of the (l+l)-dimensional NLD equation (12.1 p with Lj given 
in (12.51) in the next subsection. 

2.2. Standing wave solution 

This subsection will derive the standing wave solutions of the (1+1)- 
dimensional NLD equation (12.11) with the self-interaction (12. 5 p in the spirit 
of the technique used in 12|, [l3|, [l5|, [l7|. The solution = {ipi,ip2)'^ of the 
NLD equation (12.10 with Lj in (12. 5p . having the form 

is wanted, where m > u > 0, and |v5(x)| and |x(3^)| are assumed to decay 
exponentially as |x| — t- +oo or have the finite compact support. For such 
solution, it is not difficult to verify that both the Lagrangian and the 
energy-momentum tensor T^^ are independent of the time t, because 

Too = -i^7Vx. + m^i^ - Toi = i-07Vx, (2-9) 

Tio = -u^'j^ip, Tu = -I^'tVx + L[^], 

where i/'(x) = (^(p{x),x{x))'^ . Using the conservation laws (12.71) further gives 

Tio = -w-^tV = 0, Til = -i-^T^x + ^[*] = 0. (2.10) 
The first equation implies that ip*x is imaginary because 

-07^ = V*X + VX* = 0. 

Thus, without loss of generality, we may assume that (^{x) is real and x{^) 
is imaginary, and they are in the form 

^(^) = ( ^(^)^ = ^(^) ( cos (^(^)) \ ^ (2.11) 

where both R{x) and ^(x) are pending real functions, and R{x) is assumed 
to satisfy the inequality Q{t) = R^{x)(ix < +oo. On the other hand, 
combining the first equation in (12. 9p with the second equation in (I2.10p yields 

w-^V-w^-^^ + ^il*] = 0, (2.12) 
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which becomes for fl2.1ip 

R'^{uj -mcos{2e)) +i?2^+2^scos^+^(2^) + =0. (2.13) 
Combining (1212]) with I^M) leads to 

kujtp'''ip — kmtjjip — iipj^tp^ = 0, 
which reduces to for (12.111) 

d^ 

_ = _A;(a;-mcos(2^)). (2.14) 



Because 



I 



: tanh"^ ^ — tan (6) , 



— koj + km cos{26) ^Jk'^iw? — u"^) \ \/k'^{'m? — uj'^) 

for 6 G ( — 0.5 cos~^(aj/m), 0.5 cos~"^(a;/m)) C (— 7r/2, 7r/2) when m > w > 0, 
the solution of fl2.14p may be derived as follows: 



e{x) = tan-^ y^ ^^^ tanh (kVui"^ -u^x^j , (2.15) 

for initial data 6{0) = and k > 0. In fact, under the previous assumption, 
one may verify 6{x) G (— 7r/4, 7r/4). If coefficients s and v in (12. 5p belong to 
the set {v > 0, s > —v} for m. > a; > 0, or {f > 0, s > —v} for m > u = 0, 
then from Eq. (12.131) one has non-trivial R{x) for the localized solution as 
follows 

= ± ( "^^":i'^i"Hr" l • (2-16) 

Hereto, the standing wave solution of the NLD equation (12. ip with ( 12. 5p 
has been derived, and will be denoted as follows 

•.-i.,f)J%!fi)^^-'Ri.)(^-^';^.], (2.17) 

y V^r(x,t) J \ ism{9{x)} J 

where 6{x) and R{x) are given in Eqs. (I2.15P and (I2.16p . respectively. This 
solution represents a solitary wave with zero velocity and contains some spe- 



cial cases in the literature e.g. (16|, |47 
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Remark 2.1. It has been pointed out in |49[ that the profile of the charge 
density Jo{x,t) for the standing wave 02.171) under the scalar self-interaction 
{i.e. s ^ and f = 0) with k = 1 can be either one-humped or two-humped, 



which is also recently confirmed for any A; > by other researchers in [16 
They further pointed out there that the profile can only be one-humped for 
any A; > in the case of s = and v ^ 0. For the linear combined self- 
interaction fl2.5p with s ^ and v ^ 0, we find that the profile can also 
be one-humped (see Figs. [2] and [6] where the charge density is denoted by 
Pq(x, t) for convenience) or two-humped (see Fig. [8]). 

2.3. Solitary wave solution with nonzero velocity 

This subsection will derive another solution of the (1-1-1) -dimensional 
NLD equation (12. ip with the self-interaction (12. 5 p by using the Lorentz co- 
variance of the NLD equation. Consider a frame F with an observer O and 
coordinates {x,t). The observer O describes a particle by the wavefunction 
^{x,t) which obeys the NLD equation (12.11) with Lj given in (12. 5p . i.e. 

(^i7°^ + ^^'^- ^i^^ t) + (x, t) = 0. (2.18) 

In another inertial frame F' with an observer O' and coordinates (x', t') given 
by the Lorentz transformation with a translation in the x-direction 



X' 



' = ^{{x-Xo)-Vt), t' = ^{t- V{x - xo)), (2.19) 



which is called "boost" in the x-direction, where xq is any given position, V 
is the relative velocity between frames in the x-direction, and 7 = 1/ a/1 — 
is the Lorentz factor. According to the relativity principle, the observer O' 
describes the same particle by ^'{x',t') which should also satisfy 

(^i7°|7 + iT^^ - ^'{x', t') + {dL[[^']/ d^') (x', t') = 0. (2.20) 

Using some algebraic manipulations, the "transformation" matrix S may be 
found as follows 



^ sign(^),/^ 



2 

sign(l^)JV a/^ 



(2.21) 
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which takes ^{x,t) to under the Lorentz transformation fl2.19p . i.e. 

^'{x',t') = S^{x,t), (2.22) 

where sign(x) is the sign function which returns 1 if x > 0, if x = 0, and 
-1 if X < 0. 

Applying the transformation f l2.22p with f l2.2ip to the standing wave so- 
lution fl2.17p gives another solution of the NLD equation fl2.ip -f l23|l . i.e. the 
moving wave solution 

*™(x-Xo,t) = V 2 -/^l ^ ^ ^ 1^ 2J2 V ' M ^ ^2.23) 

This solution represents a NLD solitary wave with velocity V and will return 
to the standing wave fl2.17p if setting V = Q and Xq = 0. 

2.4. Time reversibility 

This subsection will show that the NLD equation (12. ip with Li given in 
(12. 5 p remains invariant under the time reversal operation 

x' = x, t' = -t, ^'{x',t') = K^{x,t), K = j^C, (2.24) 

where C denotes the complex conjugate operation on ^(x, t), i.e. **(x, t) = 
C^'(x,t), the time-reversal operator K satisfies 

K^K = I, K7° = 7°K, K-f^ = -j^K, (2.25) 

due to the anticommutation relation {7^,7*^} = 277^"^/, and / is the unit 
matrix. From the relations (12.250 . it can be easily verified that 

{^'^'){x',t') = (**)(x,t), 
(* V*') (x', t') = (*7°*) (x, t) , (2.26) 

(*V*')(^',0 = -(*V*)(^,t), 

so that the self-interaction Lagrangian in (12. 5p satisfies 

L[[^']ix',t') = L,[^]ix,t), (2.27) 
under the time reversal transformation (12.240 . 
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Applying the time-reversal operator K to the NLD equation fl2.18p and 
using the definition fl2.24p as well as the relations f l2.25p and f l2.27p lead 
to an equation which has the same form as shown in fl2.20p . That is, if a 
spinor ^{x, t) satisfies the NLD equation (12.180 . then the transformed spinor 
t') by the time reversal operation f l2.24p will also satisfy the same NLD 
equation, i.e. the NLD equation (12.180 is time reversible under the operation 

(EMD. 



3. Numerical methods 

As we mentioned in Section [H some numerical methods have been well 
developed for the NLD equation with a scalar or vector self-interaction. This 
section will extendedly present and compare several numerical methods for 
solving the (l + l)-dimensional NLD equation (13. ip with the general scalar 
and vector self- interaction (12. 5p . Their numerical analyses will be presented 
in Section |H 

For the sake of convenience, divide the bounded spatial domain ^2 C M 
into a uniform partition {xj\xj = jh G Q,j G Z} with a constant stepsize 
h > and give a grid in time {t„ = nr, n = 0, 1, ■ ■ ■ } with a time stepsize 
r > 0, and recast the NLD equation (12.180 into 

0^ — 

'dt ^'^'~d^ ^ ^"^^'^ ~ ^^'^'^ ~ i/v*7M*7V* = 0, (3.1) 

where 

/s := s(^+l)w^s^ M^s:=**, f.:=vik+l)wi''-'^/^, :=^Y^^7.^, 

all of which are real functions, and the dependence of ^{x,t) on {x,t) is 
implied. 

3.1. Several finite difference schemes 

Use to denote approximation of ^{xj,tn) and define the forward, 
backward and centered difference operators in space and time by: 

6^ = ±(Ef - I)/h, 5° = {E, - E-')/2h, ^ • ^ 

where / is the identity operator, and Et and E^ are the translation operators 
in time and space, respectively, defined by 



j+i' 
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whose inverses exist and are defined by 



n-1 



Besides, several symbols are also introduced for arithmetic averaging opera- 
tors: 



=1 + *"), 

and for an extrapolation operator: 



nil 

2 



4 



(70 



2 



*^ 2 J 2^ ■ 

• Crank-Nicolson schemes The CN scheme and its linearized version 
have been studied in 48|, |53|, |54] for the NLD equation with the quadric 
scalar self-interaction. For the system f l3.ip . the extension of the CN scheme 



m 



54i |62| may be written as 



(3.3) 



0, 



by approximating ( 13. ip at point with compact central difference 

quotient in place of the partial derivative, where 



r-Wi 

Fj{wj) := / fi{x)dx, I 
Jo 



S, V. 



The above CN scheme (named as CN hereafter) is fully implicit and forms a 
nonlinear algebraic system. In practice, the linearization technique is often 
used to overcome difficulty in dire ctly solving such nonlinear algebraic sys- 
tem. Two linearization techniques 62| for numerical methods of the nonlinear 
Schrodinger equation are borrowed here. The first linearized CN scheme we 
consider is using wholly the extrapolation technique to the nonlinear self- 
interaction terms in (13. ip 



St^] + 7%'A+5°*" + im^'it'f] - i^? ( /s7"* + /v*7^*7"7"* 



= 0, 
(3.4) 



12 



which will be called by LCNl. The second linearized CN scheme, denoted by 
LCN2, is 

- i£?(/.*7;.*),"7°7'^A^*," = 0, 

by partially applying the extrapolation technique to the nonlinear self-interaction 
terms. It is worth noting that the above linearized CN schemes are not lin- 
earized version of the CN scheme (13. 3p . The LCN2 scheme may conserve the 
charge and behaves better than the LCNl scheme {vide post). 

Remark 3.1. For the (l+l)-dimensional NLD equation (12. ip with the quadric 
scalar self-interaction Lagrangian (12. 3p . the CN scheme (named by CNO) pro- 



posed in j48 



IS 

St^] + j'^-fHiS^,^] + im7°^+*;^ - 2is(£+*^"^+*^")7°£+*^" = 0, (3.6) 
and its linearized version (called by LCNO) is given in [sij as follows 



2is ((**)" - rRe(*^V7'^°*")) 7°^t+*" = 0. 



(3.7) 



We will show in Section H] that the CN, CNO and LCNO schemes conserve the 
charge and the CN scheme further conserves the energy. 

• Odd-even hopscotch scheme The odd-even hopscotch scheme is a 
numerical integration technique for time-dependent partial differential equa- 
tions, see 63|. Its key point is that the forward Euler-central difference 
scheme is used for the odd grid points while at the even points the backward 
Euler-central difference scheme is recovered. Thus the odd-even hopscotch 
scheme may be explicitly implemented. Such scheme applied to the system 
(13. ip becomes 

S+^] + 7°7'5°*}' + im7°*^" 

- i£l (/,7°* + /v*7m*7°7'^*) " = 0, n + J is odd, ^^'^^ 



s^^]+' + 7^7 + i^7'*r 

n / n - n \ (3-9) 

- iil (^/s7°* + /v*7m*7%^*) . =0, n + j is even. 
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In the following we will call it by HS. 

• Leapfrog scheme The leapfrog scheme looks quite similar to the for- 
ward scheme, see e.g. fl3.8p . except it uses the values from the previous time- 
step instead of the current one. For the system f l3.ip . it is 

6^^] + 7°7'5°* ■ + im7°*," - i (/s7°* + /v*7;.*7V*)J = 0, (3.10) 

which is a three-level explicit scheme in time with a central difference in 
space and will be named by LF. 

• Semi-implicit scheme Another three-level scheme considered here for 
the system (13.11) is 

S^^]Wl'iX^] + im-f'i'i^^-i (/,7°* + /v*7^*7V*)" = 0, (3.11) 

which is obtained by approximating explicitly the nonlinear terms but im- 
plicitly the linear terms and will be called by SI. It is worth noting that 
such semi-implicit scheme has been studied for the NLD equation with the 
quadric scalar self-interaction in 46 



3.2. Exponential operator splitting scheme 

This subsection goes into discussing exponential operator splitting scheme 
for the NLD equation (13. ip . For convenience, we rewrite the system (13. ip as 
follows 

^t = iC + K+K)^, (3.12) 

where £, A/^ and AC are the linear operator, and nonlinear operators for 
self-interaction terms, respectively. That is, they are defined by 

Cqi ■= -^o^^q,^ - im7°*, A^* := i/s7°*, A/;* := i/v*7^*7V*- 

Then the problem (I3.12p may be decomposed into three subproblems as fol- 
lows 

*t = C^, (3.13) 
*i=A4*, (3.14) 
^t=K^- (3.15) 
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Based on the exact or approximate solvers of those subproblems, a more 
general i^-stage A^-th order exponential operator splitting method 6J] for 
the system fl3.12p can be cast into 



K 



= X{{^Mr^J^■'')^Mr^^^■^)^Mr^^^t^))'^^l. (3.16) 



i=l 



where Tj denotes the time stepsize used within the i-th stage and satisfies 
Yl!i=i'Ti = ''"5 ^'^d {^j-"^'', ^j-^'', ^P^} is any permutation of {£, A/'s, A/'v}. Here- 
after we call the operator splitting scheme fl3.16p by OS(A^). 

A simple example is the well-known second-order accurate operator split- 
ting method of Strang (named by OS (2)) with 

K = 2, Ti = T2, A?=A^\ Af=A^\ Af'^=A^\ (3.17) 

Another example is the fourth-order accurate operator splitting method [g 
with 

i^=8, A^^^=Af\ Af=A^^\ Af=Af\ g = 1,4,6, 7,p = 2,3,5, 



r 7+ V13- a/2(1 + V13) 

ri =T8 = , T2 = T7= — r, 

5- VT3 + ^J2{l + ^/T3) 

T3 =T6 = , T-4 = Ts 



T2 - Ti 3ti - 2T2 

which is denoted by OS (4) in the following. 



Remark 3.2. Another operator splitting scheme is studied in [52| for the 
NLD equation (12. ip but only with the quadric scalar self-interaction La- 
grangian, and the second-order accurate Strang method (I3.17P is applied 
there. For the system (13. ip . it is based on the following operator decompo- 
sition 

(£+A4+k)*, (3.18) 

with 

:= -7V*x, >^s* := "i {m - /s) 7°*, K"^ ■= i/v*7;.*7V*- 
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3.2.1. Linear subproblem 

We are now solving the the hnear subproblem f l3.13p . Denote its "initial 

data" by = {{^pi)f\ Mff at the z-th stage in i^M)- 

If the spinor ^ is periodic {e.g. 27r-periodic) with respect to x, the Fourier 
spectral method is employed to solve fl3.13l) and gives 

= ^-i(^exp (-ir,(/s:7V + ^7°))^ (^f) )• (3.19) 

Here J-" and J-""^ denote the discrete Fourier transform operator and its in- 
verse, respectively, defined by 



(^(*))^ := *i exp (-i27rK^) k = 0, . . . , J - 1, 

1=0 ^ 

(^-i(*))^. :=- J2 *K exp (i27rj^) j = 0, . . . , J - 1, 

where J is the grid point number, and the matrix exponential in f l3.19p can 
be easily evaluated as follows 



exp (— irj(K7*^7^ + 7777°)) 



cos(Cri) - sin(Cri) -i| sin(Cri) 

-if sin(Cri) cos(Cri) + sin(Cri; 

(3.20) 
with ( = a/k^ + m?. 

When the spinor is not periodic with respect to x, the fifth-order 
accurate finite difference WENO scheme will be used to solve the linear 
subproblem fl3.13p . The readers are referred to 66|] for details. In this case, 
the linear subproblem (13.130 can also be solved by using the characteristics 
method. 



3.2.2. Nonlinear subproblems 

The nonlinear subproblems f l3.14p and f l3.15p are left to be solved now. 
Their "initial data" is still denoted by = {{4^i)f\ (^2)^^)^ at the z-th 
stage in fl3.16p . and define 

i-l 

^li^ tn ~^ ^ ] Tp; i 1; 2, ■ ■ ■ , K. 
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For nonlinear subproblem fl3.14p . it is not difficulty to verify that dfWs = 
so that 

dtfs = 0. (3.21) 

Using this local conservation law gives the solution at t = of (jSH]) 

with the "initial data" as follows 

,(0) 



*r = 1^ (/s).7°dt j *f = exp (i(/s)f 7°r.) 

= diagjexp (i(/s)f r,) ,exp (-i(/3)f r,) } *f . (3.22) 

For the nonlinear subproblem fl3.15p . one may still similarly derive the 
following local conservation laws 

dti^Jo^) = 0, 5t(*7i*) = 0, dtU = 0. (3.23) 

by direct algebraic manipulations if using the fact that ^-/q^, "^-f^^ and 
/v are all real. Consequently, integrating fl3.15p from til^ to til^^'' gives its 
solution as follows 

= exp (i(A*7^*)f tS'^t.) 

(■ \( cos(/3r,) isin(/3ri) \ (0) 
= exp(iar,) \ . . ,^ . * .^ (3.24) 

Y ism(/3ri) cos(/3ri) j 

where a = (/v^To*)^ ^"^^ = (/v*7i*)f • 

Remark 3.3. It is because the local conservation laws fl3.2ip and (I3.23P 
are fully exploited here that we can solve exactly the nonlinear subproblems 
fl3.14p and (13.150 which imply the more higher accuracy of the OS method 
than that of other methods. 

In summary, we have 

• The CN (13. 3p and CNO (13. 6 P schemes are nonlinear and implicit, and 
could be solved by iterative algorithms such as Picard iteration and 
Newton method. 

• The LCNO (ETD, LCNl LCN2 ((331) and SI (131111 schemes are linear 
and implicit. 

• The HS dSlD-dSlD, LF ^AO\f . and OS(A^) fl3l6l) schemes are explicit. 
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4. Numerical analysis 



Before investigating the performance of the numerical methods proposed 
in Section [3l this section will go first into numerical analysis of them, includ- 
ing the accuracy in the sense of the truncation error, time reversibility and 
the conservation of the charge or energy. 

Proposition 4.1. G C°°(M x [0,+cx))) is periodic, then the CN, 

CNO, LCNO, LCNl, LCN2, HS, LF and SI schemes are of order 0{t'^ + h^), and 
the OS(A^) scheme is of order 0{t^ + h"^) for any arbitrary large m > 0. 

Proof. The proof is very straightforward by using directly the Taylor series 
expansion for the finite difference schemes and the Fourier spectral analysis 
for the OS(A^) scheme, and thus is skipped here for saving space. □ 

Proposition 4.2. The CN, CNO, HS, LF, SI, and OS(iV) schemes are time 
reversible, but the LCNO, LCNl, LCN2 schemes are not. 

Proof. We give the proof for the CN and LCNl schemes as an example and 
the others can be proved in a similar way. 

According to the transformation fl2.24p . the relation between the trans- 
formed finite difference solution and the original one should be — 
= (i^^')" with n' = —n. Consequently, we have 

K6+^" = i- = = -6ti^T '\ 

r r (4 1) 

t ] 2 2 ^ ^ '3 



and then using the relations in fl2.26p yields 

-(*7i*)," = -^+(*'7i*') • 



+ (*7on" = ^*^(*7o*')-'-\ 

l\n'-l 



for / G {s, v}. Applying the time- reversal operator K to the CN scheme (13.31) 
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and using the commutation relation fl2.25p and Eqs. fl4.1l) and fl4.2p lead to 



which is exactly the CN scheme fl3.3p applied to That is, the CN 

scheme is invariant under the the time-reversal transformation, namely, it is 
time reversible. 

The fact that the LCNl scheme (13 ■4p is not time reversible can be observed 
directly if noting 

- Kit, (/s7°* + /v*7p*7°7"*) 

3 / ] / \ n'+i 

I (/sV*' + /;*'7,* 7V*')^. - 2 (/sV*' + /;*'7,*V7^*';^ 

n'-l 



1 



^it, (/3'7°*' + />'7m*'7°7^*' 



□ 



Next, we will discuss the conservation of the discrete energy, linear mo- 
mentum and charge defined below for the numerical methods given in Section 
[3l After performing the integration in the computational domain Q = [xl, xr] 
and then approximating the first derivative operator dx with the centered 
difference operator as well as the integral operator J^^ dx with the sum- 



mation operator /iX]j=i -'^l- P-Sp . we have the discrete energy, linear 



momentum and charge at the n-th time step 
J 



J 

P," = /.^Im(*t50^)^n 
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where the inner product (■, ■) is defined as 



for two complex-valued vectors u and v, and J is the grid point number. 
Here the values of * at Xj with 1 < j < J are unknowns and those at Xq 
and xj+i are determined by appropriate boundary conditions. 

We first present the following lemma which can be verified through direct 
algebraic manipulations and will be used in discussing the conservation of 
the discrete charge, energy and linear momentum. 

Lemma 4.3. Given that is a complex-valued vector mesh function with 
two components evaluated at the mesh {xj,tn} (n = 0,1, ■ ■ ■ , N , j = 0,1, ■ ■ ■ , J+ 
1) and a matrix T G {/2,7°,7^j ^l^J^}, we have the following identities 



(d) 2Re((5+*^"r£+*^") = 5+(*^r*j"); 

(e) 2ilm (7°7H£+5S)*", 5+*") = -5+ {h E/=i(*^''5S*,")) 



holds for any two spinors u{x, t),v{x, t) and a G {+, — , 0}, and then we have 




Proof. It can be checked that the following Leibniz rules 




J J 




Thus the identity (a) holds. 
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Because the operator —S^ is the adjoint operator of and 7^7^ is an 
Hermite matrix, we get the identity (6) directly by rearranging the summa- 
tion. The identity (c) can be easily verified if using the fact (7°r)''' = ■y^T. 
The proof of (d) (resp. (e)) is similar with that of (6) (resp. (c)). □ 

Proposition 4.4. The CN, CNO, LCNO, LCN2, and OS (TV) schemes conserve 
the discrete charge, but only the CN scheme conserves the discrete energy. 

Proof. We begin with the discrete conservation law of charge for the CN 
scheme ( 13. 3p . Performing the inner product of if^"" and the CN scheme (13. 3p 
leads to 



..0/ 



(4.3) 



i.^^A+(*7;.*)"7V^t+*" 



0, 



and then the conservation law of the discrete charge can be easily verified by 
taking directly the real part as follows 

(4.4) 

where Lemma|13](a) is applied to the first term in Eq. 04.31) . (b) to the second 
term and (c) to the third and fourth terms. Similarly, it can be verified that 
fOD holds for the CNO ([M]), LCNO ([STD and LCN2 schemes. 

Performing the inner product of the CN scheme (13.31) and 5^**^, keeping 
the imaginary part and applying Lemma 14.31 give directly the conservation 
of the discrete energy 



^tE'f, = -Im [{n^'}^,i'5t^-j + £+*}V5+*}+i) 



For the Fourier spectral method (I3.19p . using the fact that exp (— irj(K7°7^ + m7°)) 
in (I3.20p is a unitary matrix yields 



exp (-irj(h;7°7^ + m7°)) 



(1) 



2 






3 
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and then 



II' = ll^-i [exp (-ir,(K7V + "^7°)) ^ if 

= ||exp(-ir,(«:7V + m7°))-F(*W)|f = ||^(*W)||' = ||*(°)|f , 

where Parseval's identity is apphed twice. It can be readily verified that 
the matrix exponents in Eqs. fl3.22p and fl3.24|] are unitary, thus H^'^^-'p = 
holds both for Eqs. ( K^ and (^^, i.e. Qt should be conserved 
for solutions of the nonlinear subproblems. Therefore, the OS(A^) scheme 
satisfies the conservation law of charge. □ 

Remark 4.1. It will be verified by numerical results in Section |5] that the 
LCNl, HS, LF and SI schemes do not conserve the discrete charge or energy, 
and none of the numerical methods presented in Section [3] conserves the 
discrete linear momentum. 



5. Numerical results 

This section will conduct numerical simulations to compare the perfor- 
mance of numerical schemes proposed in Section [3] and then utilize the OS (4) 
scheme to investigate the interaction dynamics for the NLD solitary waves 
fl2.23p under the scalar and vector self-interaction. For those localized NLD 
solitary waves, the periodic boundary condition for the OS(iV) scheme and 
the non-refiection boundary condition for other schemes could be adopted at 
the boundaries of the computational domain if a relatively large computa- 
tional domain has been taken in our numerical experiments. 

All calculations are performed on a Lenovo desktop computer with Intel 
Core 15 650 CPU and 4GB RAM using double precision in the 3.0.0-24- 
generic x86_64 Linux operation system and the compiler is gcc 4.6.1. The 
computational domain Q will be taken as [—50,50] in Examples 15. ItlSTSl and 
[—100, 100] in Example 15.61 and the particle mass m in Eq. fl2.18p is chosen 
to be 1. 

Example 5.1. The first example is devoted to comparing the numerical per- 
formance of all the numerical methods in Section [3] in terms of the accuracy, 
the conservativeness, the efficiency and the error growth. A one-humped 
solitary wave with the velocity V = —0.2 is simulated here under the quadric 
scalar self-interaction [i.e. v = and k = 1), traveling from right to left with 
the parameters in ([22S]): xq = 5, s = 0.5, and u = 0.75. The P^-RKDG 
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method [50[ is also included here for comparison, which is assembled with a 
fourth-order accurate Runge-Kutta time discretization in time and the Leg- 
endre polynomials of degree at most N as local basis functions in the spatial 
Galerkin approximation. 

Tables [1] and [2] summarize the numerical results at the final time t = 50, 
where err2 and erroo are respectively the and l°° errors at the final time, 
Vq, Ve, Vp measure respectively the variation of charge, energy and lin- 
ear momentum at the final time relative to the initial quantities, and the 
CPU time of calculations with the same mesh size is recorded for comparing 
the efficiency. It can be observed clearly there that: (1) while the CN, CNO, 
LCNO, LCNl, LCN2, HS, LF, SI and OS (2) schemes are of the second-order 
accuracy, the OS (4), P^-RKDG and P^-RKDG methods exhibit at least 
the fourth-order accuracy; (2) The CN, CNO, LCNO, LCN2, OS (2) and OS (4) 
schemes conserve the discrete charge and only the CN scheme conserves the 
discrete energy, but none conserves the discrete linear momentum; (3) The 
OS (4) scheme could also keep very accurately the discrete energy and linear 
momentum with relatively fine meshes. All above numerical results are con- 
sistent with the theoretical results given in Section HI Among the numerical 
methods of the second order accuracy, it is also found that the OS (2) scheme 
runs fastest (8.2 seconds for the mesh r = 0.005 and h = 0.20) if requiring 
to attain almost the same accuracy. Similarly, the OS (4) scheme runs much 
more faster than both P^-RKDG and P^-RKDG methods, and the ratio of 
the CPU time used by the OS (4) scheme over that used by the P^-RKDG 
method is around 3.05%, and reduces to around 2.26% over that used by the 
P^-RKDG method. 

Fig. [1] plots the error history in the finest mesh used in Tables [1] and 
[2J According to the curves shown there, it can be seen there that the l°° 
error of all the schemes increases almost linearly with the time. However, 
the slopes, obtained by the linear fitting, are different. The smaller the slope 
is, the longer time the scheme could simulate to. The SI scheme has the 
largest slope 1.412 x 10^*^'^ while the OS (2) scheme has the smallest one 
2.334 X 10~°^ among all the second-order accurate methods. Further, the 
slopes of the curves of errors for the OS (4), P^-RKDG and P'^-RKDG 
schemes are almost the same value of 3.199 x 10"^'^ which is much more 
smaller than those of the second-order accurate schemes. 

Remark 5.1. Both the theoretical and numerical comparison of the OS (2) 
scheme with the CNO and LCNO schemes show that the former is better, espe- 
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Table 1: Example 15. II Part I: Numerical comparison of the accuracy, the conservativeness 
and the efficiency at t = 50 with the time stepsize being set to r = ^h. The CPU time in 
seconds is recorded for the finest mesh. 





h 


Vq 


Ve 


Vp 


err2 


Order 


err^ 


Order 


Time (s) 


CN 


0.08 


1.38E-1.5 


1.09E-15 


2.6.5E-06 


2.74E-02 




2.61E-02 








0.04 


4.66E-1.5 


1.37E-15 


1.6.5E-07 


6.84E-03 


2.00 


6.50E-03 


2.00 






0.02 


4.78E-15 


1.37E-1.5 


1.03E-08 


1.71E-03 


2.00 


1.63E-03 


2.00 






0.01 


3.05E-14 


2.05E-15 


6.42E-10 


4.27E-04 


2.00 


4.06E-04 


2.00 


367.1 


CNO 


0.08 


7.0.5E-15 


3.24E-08 


2.37E-06 


2.31E-02 




2.21E-02 








0.04 


l.OlE-1.5 


2.02E-09 


1.47E-07 


5.76E-03 


2.00 


5.51E-03 


2.00 






0.02 


5..54E-1.5 


1.26E-10 


9.17E-09 


1.44E^03 


2.00 


1.38E-03 


2.00 






0.01 


2.88E-14 


7.88E-12 


5.72E-10 


3.59E-04 


2.00 


3.44E-04 


2.00 


345.3 


LCNO 


0.08 


2.14E-15 


1.83E-06 


4.31E-05 


2.75E-02 




2.62E-02 








0.04 


1.26E-16 


2.29E-07 


5..55E-06 


6.87E-03 


2.00 


6.53E-03 


2.00 






0.02 


8.81E-15 


2.86E-08 


7.04E-07 


1.72E-03 


2.00 


1.63E-03 


2.00 






0.01 


3.11E-14 


3.58E-09 


8.87E-08 


4.29E-04 


2.00 


4.08E-04 


2.00 


81.8 


LCNl 


0.08 


2.22E-04 


1.92E-04 


3.73E-04 


3.53E-02 




3.34E-02 








0.04 


2.78E-0.5 


2.40E-05 


4.6.5E-05 


9.06E-03 


1.96 


8.56E-03 


1.97 






0.02 


3.47E-06 


3.01E-06 


5.80E-06 


2.30E-03 


1.98 


2.17E-03 


1.98 






0.01 


4.34E-07 


3.76E-07 


7.2.5E-07 


5.78E-04 


1.99 


5.45E-04 


1.99 


118.7 



LCN2 


0.08 


2.77E-15 


1.64E-07 


7.01E-06 


2.77E-02 




2.64E-02 








0.04 


7..55E-16 


1.98E-08 


6.77E-07 


6.92E-03 


2.00 


6.58E-03 


2.00 






0.02 


7.55E-16 


2.44E-09 


7.24E-08 


1.73E-03 


2.00 


1.64E-03 


2.00 






0.01 


2.91E-14 


3.03E-10 


8.29E-09 


4.32E-04 


2.00 


4.11E-04 


2.00 


118.4 


HS 


0.08 


1.62E-06 


1..55E-06 


4.33E-06 


2.00E-02 




1.5.5E-02 








0.04 


l.OlE-07 


9.65E-08 


2.70E-07 


.5.00E-03 


2.00 


3.88E-03 


2.00 






0.02 


6.30E-09 


6.03E-09 


1.68E-08 


1.25E-03 


2.00 


9.71E-04 


2.00 






0.01 


3.94E-10 


3.77E-10 


1.0.5E-09 


3.13E^04 


2.00 


2.43E-04 


2.00 


14.8 


LF 


0.08 


5.88E-05 


4.73E-05 


8.68E-06 


1.41E^02 




1.35E-02 








0.04 


5.51E-06 


4.43E-06 


6..59E-07 


3.51E-03 


2.01 


3.36E-03 


2.01 






0.02 


6.24E-07 


5.01E-07 


6.91E-08 


8.74E-04 


2.00 


8.36E-04 


2.00 






0.01 


7.54E-08 


6.05E-08 


8.19E-09 


2.18E-04 


2.00 


2.09E-04 


2.00 


8.9 


SI 


0.08 


3.79E-08 


1.15E-07 


3.94E-06 


3.59E-02 




4.76E-02 








0.04 


3.68E-09 


8.36E-09 


2.40E-07 


8.97E-03 


2.00 


1.19E-02 


2.00 






0.02 


4.09E-10 


6.87E-10 


1.44E-08 


2.24E-03 


2.00 


2.98E-03 


2.00 






0.01 


4.89E-11 


6.46E-11 


8.20E-10 


.5.60E-04 


2.00 


7.44E-04 


2.00 


132.7 
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Tabic 2: Example 1 5. II Part II: Numerical comparison of the accuracy, the conservativeness 
and the efficiency at t = 50. The CPU time measured in seconds is listed for the finest 
mesh. 



T h Vq Ve Vp err2 Order err^ Order Time (s) 

OS (2) 0.04 0.78 1.04E-13 2.33E-07 3.97E-05 2.03E-03 9.96E-04 

0.02 0.39 1.33E-13 1.45E-08 2.55E-13 2.53E-04 3.01 1.90E-04 2.39 

0.01 0.39 3.58E-13 9.04E-10 6.08E-13 6.33E-0.5 2.00 4.7.5E-05 2.00 

0.005 0.20 6.15E-13 5.60E-11 7.13E-13 1.58E-05 2.00 1.21E-05 1.97 8.2 

OS (4) 0.04 0.78 7.40E-13 5.78E-13 4.12E-0.5 1.71E-03 4.27E-04 

0.02 0.39 1.46E-12 1.25E-12 2.17E-12 7.92E-08 14.40 3.32E-08 13.6-5 

0.01 0.20 1.88E-12 1.60E-12 2.71E-12 4.28E-10 7.53 3.28E-10 6.66 

0.005 0.10 7.00E-13 6.55E-13 2.44E-12 1.93E-11 4.47 1.44E-11 4.51 16.8 

P3_RKDG 0.04 0.78 1.03E-05 5.47E-07 2.61E-06 1.46E-04 1.44E-04 

0.02 0.39 8.60E-08 5.96E-07 6.31E-07 6.59E-06 4.47 8.85E-06 4.03 

0.01 0.20 6.98E-10 4.32E-08 4.69E-08 4.10E-07 4.01 5.62E-07 3.98 

0.005 0.10 7.88E-12 2.75E-09 3.00E-09 2.56E-08 4.00 3.69E-08 3.93 551.6 

p-^-RKDG 0.04 0.78 l.OOE-07 2.13E-07 l.lOE-07 8.02E-06 1.04E-05 

0.02 0.39 8.53E-10 4.51E-09 3.44E-09 2.58E-07 4.96 3.36E-07 4.96 

0.01 0.20 2.28E-11 6.41E-11 5.22E-11 8.46E-09 4.93 1.16E-08 4.86 

0.005 0.10 2..58E-12 9.56E-11 2.16E-10 3.09E-10 4.78 2.39E-10 5.60 744.8 



cially in terms of efficiency and error growth. Therefore in some sense this 
is an answer to the debate stimulated in 52|, |53| over twenty years ago on 



which one is most efficient among the OS (2), CNO and LCNO schemes. 

Example 5.2. The P^-RKDG method has been successfully applied be- 
fore into investigating the interaction for the NLD solitary waves under the 



quadric scalar self-interaction in |49|-|51| , but the numerical comparison shown 



in Example 15. II tells us that the proposed OS (4) scheme should be preferred 
now. In this example, we further conduct numerical comparison among the 
OS (4), P^-RKDG and P^-RKDG methods in simulating one-humped and 
two-humped solitary waves. Two typical profiles of the charge density for 
the NLD solitary wave displayed in Fig. [2] are considered, one denoted by 
Case 1 has a two-humped proffie under the quadric scalar self-interaction, 
and the other denoted by Case 2 has a one-humped profile under the cubic 
scalar and vector self-interaction. These two solitary waves are located ini- 
tially at Xq = 5, travel from right to left with the velocity V = —0.2 and stop 
at the final time t = 50. For convenience, we use pq^x.t) = Jq to represent 
the charge density. 

The numerical comparison is shown in Table El from which we can observe 
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f 

Figure 1: Exainple l5.ll The l°° error history. 

Table 3: Example [5?2] Numerical comparison among the OS (4) , P^-RKDG and P^-RKDG 
methods. The CPU time is measured in seconds. 

Case 1 in the mesh of r = 0.01 and h = ^ 







Ve 


Vp 


eiTa 


eri-oo 


Time (s) 


OS (4) 


1.96E-12 


l.lOE-12 


3.45E-12 


2.15E-09 


2.17E-09 


7.8 


P-^-RKDG 


3.34E-08 


2.27E-07 


2.96E-07 


2.89E-06 


4.15E-06 


146.6 


P^-RKDG 


1.35E-12 


9.43E-08 


5.34E-08 


9.02E-08 


6.94E-08 


195.5 


Case 2 in the mesh of r = 


0.005 and h = 


100 
1024 




Vq 


Ve 


Vp 


eiTa 




Time (s) 


OS (4) 


9.82E-13 


8.35E-13 


2.95E-12 


7.33E-11 


8.63E-11 


46.8 


P-'-RKDG 


3.21E-10 


9.04E-09 


1.95E-08 


2.00E-07 


3.98E-07 


881.1 


P^-RKDG 


9.10E-13 


4.89E-09 


1.57E-09 


4.68E-09 


6.35E-09 


1156.5 



that, (1) with the same mesh, no matter sparse (r = 0.01 and h = ^ ^ 
0.1953) or fine (r = 0.005 and /i = ^ ^ 0.0977), the OS (4) scheme is 
more conservative and higher accurate than both P^-RKDG and P^-RKDG 
methods; (2) the OS (4) scheme runs much faster than both P^-RKDG and 
P^-RKDG methods as we have found in Table [2l Here, the ratio of the CPU 
time used by the OS (4) scheme over that used by the P^-RKDG method is 
around 5.32%, and reduces to around 4.05% over that used by the P^'-RKDG 
method for both cases. The error history is plotted in Fig. |3l which shows 
that those three methods all have almost zero slope (under 4.50E-09). This 
can be used to explain our previous success of the P'^-RKDG method in (igj - 



5l| . Further more, the more smaller errors of the OS (4) method mean that 
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Figure 2: Exainple l5.2l The initial charge density pq{x, t) for two typical cases, V = —0.2, 
xq = 5. Case 1 is shown in the solid line, a two-humped profile (w — 0.3) under the quadric 
scalar self-interaction (fc = 1, s = 0.5 and v = 0); Case 2 is shown in the dashed line, a 
one-humped profile {uj ~ 0.75) under the cubic scalar and vector self-interaction (fc = 2 
and s = V = 0.5) 

it should be more powerful than others. 



to"' 




Case 1 




10"" 


4.322 X 10^ 


DS(4) 

^P-'-RKDG 


8 -. 

^ 10 






5.012 X 10"" 




10' 


3.89x10"" 


10"^° 






5 


10 15 20 25 30 35 
/ 

Case 2 


40 45 5 


10"^ 


5.376 X 10-'° 


03(4) 

^^P^-RKDG 
^P*-RKDG 








-2.653 X 10"" 




I0"^° 


1.533 X 10-" 







5 10 15 20 25 30 35 40 45 50 

f 



Figure 3: Example 15.21 The l°° error history. The slopes are displayed above the curves. 

It has been shown that the OS (4) scheme behaves best for both one- 
humped and two-humped NLD solitary waves in long time simulations. There- 
fore, we conclude the comparison with the judgement that the OS (4) scheme 
is the most suitable for simulating the interaction dynamics for the NLD 
solitary waves in terms of the accuracy, the conservativeness, the efficiency 
and the error growth. The OS (4) scheme will be utilized to investigate the 
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binary collision of the NLD solitary waves. The initial setup is the linear 
superposition of two moving waves 

qf{x,t = 0) = *[""(x-xi,t = 0) + *[""(x-Xr,t = 0), (5.1) 

where *™™(a; — Xpos, t) denote the moving waves f l2.23p centered at Xpos with 
the speed Vpos and the frequency Upos for pos G {1, r}. In the following 
examples, two equal solitary waves are placed symmetrically at t = with 
—xi = = 10 and V\ = — VJ. = 0.2. Several typical NLD solitary waves 
are considered with the parameters given in Table HI and both quadric {k = 
1) and cubic {k = 2) cases will be studied. It should be noted that the 
experiments carried out in the literatures are all limited to the collisions of 
the NLD solitary waves under the quadric scalar self-interaction. A relatively 
fine mesh, r = 0.005 and h = 100/2^^ ^ 0.0122, is adopted hereafter. 



Table 4: The initial setups of different cases in binary collisions. 



case 


s 


V 




Remarks 


Bl 


0.5 





0.8 


scalar, one- humped 


B2 





0.5 


0.8 


vector, one-humped 


B3 


0.5 


0.5 


0.8 


scalar and vector, one-humped 


B4 


0.5 





0.3 


scalar, two-humped 


B5 


4 


1 


0.1 


scalar and vector, two-humped 



Example 5.3. The collision of two equal one-humped solitary waves under 
the scalar self-interaction, i.e. Case Bl in Table HI is studied in this example. 
The interaction dynamics for the quadric case are shown in the left plot of 
Fig. m where two equal waves with the initial amplitude of 0.4082 move 
close at a velocity of 0.2 and overlap each other, then separate into a left 
moving wave and a right moving wave with the amplitude of 0.3743 and the 
velocity of 0.1831. Similar phenomena are observed for the cubic case shown 
in the right plot of Fig. H] except that (1) two waves overlap more stronger 
around t = 41 now due to the stronger nonlinearity; (2) after collision, the 
amplitude decreases to 0.5899 from the initial amplitude of 0.6455 while the 
velocity also decreases to 0.1037. In both cases, the discrete charge, energy 
and linear momentum are approximately conserved in the interaction since 
the variation of them at t = 80 is under 1.58E-10. 
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(a) fc = 1 



(b) A: = 2 



Figure 4: Example 15.31 Binary collision of the NLD solitary waves under the scalar self- 
interaction. 



Example 5.4. The collision of two equal one-humped solitary waves under 
the vector self-interaction, i.e. Case B2 in Table HI is studied in this example. 
To the best of our knowledge, it is the first time to study binary collision 
of the NLD solitary waves under the vector self-interaction. The interaction 
dynamics for the quadric case are shown in the left plot of Fig. O where the 
waves keep the shape and the velocity after the collision. A totally different 
phenomenon appears for the cubic vector self-interaction as displayed in the 
right plot of Fig. [51 The initial one-humped equal waves first merge into 
a single wave, then separate and overlap again. Around t = 50, collapse 
happens and highly oscillatory waves are generated and moving outside with 
a big velocity near 1, meanwhile a one-humped wave with small amplitude 
is formed at the center. In both cases, the discrete charge, energy and linear 
momentum are approximately conserved in the interaction since the variation 
of them at t = 100 is under 5.4115-11. Note in passing that the collapse here 



is different from that shown in [49[. It was reported there that the strong 
negative energy and radiation appear when the collapse happens during the 
binary collision of two-humped waves. 

Example 5.5. This example is devoted into investigating for the first time 
the collision of two equal NLD solitary waves under the scalar and vector 
self-interaction, i.e. Case B3 in Table HI The interaction dynamics for the 
quadric case are shown in the left plot of Fig. [6l where two equal waves with 
the initial amplitude of 0.2041 move close at a velocity of 0.2 and overlap 
each other, then separate into a left moving wave and a right moving wave 
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Figure 5: Example 15.41 Binary collision of the NLD solitary waves under the vector self- 
interaction. 



with the amphtude of 0.2091 and the velocity of 0.1968. The collapse similar 
to that shown in right plot of Fig. happens again for the cubic vector self- 
interaction, see the right plot of Fig. [61 The initial one-humped equal waves 
first merge into a single wave at t = 38, then separate and overlap again. 
Around t = 50, collapse happens and highly oscillatory waves are generated 
and moving outside with a big velocity near 1. In both cases, the discrete 
charge, energy and linear momentum are approximately conserved in the 
interaction since the variation of them at t = 80 is under 3.53E-10. 



Example 5.6. As reported before in |49|, [SlJ, collapse happens in binary 
and ternary collisions of the NLD solitary waves under the quadric scalar 
self-interaction if the two-humped waves are evolved. In this example, we 
will show further that collapse could happen in binary collision of equal two- 
humped waves under the cubic scalar self-interaction and under the linear 
combination of scalar and vector self-interactions. First, Case B4 in Table H] 
is studied and the interaction dynamics are shown in Fig. [TJ which clearly 
shows that (1) collapse happens in both quadric and cubic cases but is more 
stronger in the latter; (2) two initial waves at the same velocity are decom- 
posed into groups with different velocities after the collision, but there is 
no such decomposition for the cubic case. In the left plot of Fig. [71 the 
highly oscillating waves with small amplitude move outside at a big velocity 
of 0.9644, while the one- humped waves with big amplitude follow them at a 
small velocity of 0.4626. In both cases, the discrete charge, energy and linear 
momentum are approximately conserved in the interaction since the variation 
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(a) fc = 1 



(b) k = 2 



Figure 6: Example 15.51 Binary collision of the NLD solitary waves under the scalar and 
vector self-interaction. 

of them at t = 80 is under l.OlE-5. Second, binary collision of equal two- 
humped solitary waves under the scalar and vector self-interaction, i. e. Case 
B5 in Table HI is plotted in Fig. [HI The phenomena are very similar to that 
shown in Fig. [TJ and the "decomposition" phenomenon for the quadric case 
is more obvious than that shown in the left plot of Fig. [71 

6. Conclusion and outlook 

Several numerical methods for solving the NLD equation with the scalar 
and vector self-interaction have been presented and compared theoretically 
and numerically. Our results have revealed that among them, the OS (4) 
scheme, one of the fourth-order accurate OS methods, performs best in terms 
of the accuracy and the efficiency. Particularly, the OS (4) scheme is usually 
more accurate than the P^-RKDG method in the mesh of the same size, but 
the former needs much more less computational cost than the latter. Such 
superior performance of the OS methods is credited to the full use of the local 
conservation laws of the NLD equation such that the nonlinear subproblems 
resulted from them are exactly solved. The interaction dynamics for the 
NLD solitary waves under the quadric and cubic self-interaction have been 
investigated with the OS (4) scheme. We have found that such interaction 
dynamics depend on the exponent power of the self-interaction. Actually, it 
has been observed for the first time in our numerical experiments that, (1) 
collapse happens in collision of two equal one-humped NLD solitary waves 
under the cubic vector self-interaction but such collapse does not appear 
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5- 8s 




(a) fc = 1 {h) k = 2 



Figure 7: Example 15.61 Binary collision of the two- humped NLD solitary waves under the 
scalar self-interaction. 

for corresponding quadric case; (2) two initial waves at the same vefocity 
are decomposed into groups with different vefocities after collapse in binary 
collision of two-humped NLD solitary waves under the quadric scalar self- 
interaction or under the quadric scalar and vector self-interaction but such 
phenomenon does not show up for corresponding cubic case. More efforts 
on the interaction dynamics for the NLD solitary waves under more general 
self-interaction with the OS (4) method are still going on. 
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